ct = data4d('../data/lung/ct', 'short')
itv = data3d('../data/lung/itv', 'short')
step = float_param(0.002, 0.001, 0.01, 'S')

# rpl params
colormap = data1d_rgba('jet.png', 'uchar')
table = data1d('density-to-rpl-table.png', 'uchar')
rpl_max = float_param(0.8, 0, 4, 'M')

# transfer function params
tf = data1d_rgba('bone.png', 'uchar')
tf_pos = float_param(0, -200, 200, 'P')
tf_width = float_param(100, 10, 400, 'W')
tf_opacity = float_param(0.8, 0, 1, 'O')

# lighting params
ambient = float_param(0.1, 0, 1, 'A')
lightpos = float3_param(-3, 3, 3)

cut = float_param(-2, -2, 2, 'C')

C = normalize(S - E)
dP = (-C) * step
ray_length = length(E - S)
jitter = rand() * step
    
# compute frame
frame = fmod(T*5, 10)
    
# first pass: find distal edge of ITV
t = jitter
distal = -1.0
P = S

while t < ray_length:
    q = linear_query_3d(itv, P)
    if q > 0: distal = t
    P += dP
    t += step

# second pass: compute RPL
P = S
t = jitter
rpl = 0.0

while t < distal:
    density = linear_query_4d(ct, P, frame) * 32768
    rel_density = 0.21 + (1.78 - 0.21) * linear_query_1d(table,
        (density - (-766)) / (1173 - (-766)))
    rpl += rel_density * step
    P += dP
    t += step
    
rpl_color = float3(1)
if rpl > 0:
    rpl_color = linear_query_1d_rgba(colormap, rpl / rpl_max).xyz

# third pass: volume render
t = jitter
P = S
result = float3(0)
see_through = 1.0

while t < ray_length:
    density = cubic_query_4d_cut(ct, P, C, cut, frame) * 32768
    tf_query = (density - tf_pos + tf_width) / (2 * tf_width)

    if tf_query >= 0:
        rgba = linear_query_1d_rgba(tf, tf_query)

        color = rgba.xyz * rpl_color
        # if rpl > 0: color = rpl_color
        
        N = -normalize(cubic_gradient_4d_cut(ct, P, C, cut, frame))
        L = normalize(lightpos - P)
        lit = phong(L, N, C, color, 0.5, 50, ambient)

        result += see_through * rgba.w * tf_opacity * lit
        see_through *= 1 - rgba.w * tf_opacity
        if see_through < 0.01: break

    P += dP
    t += step

return result
